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We consider an efficient scheme to simulate fermionic Hubbard models with nonlocal density- 
density interactions in two dimensions, based on bond-centered auxiliary-field quantum Monte Carlo. 

The simulations are shown to be sign-problem free within a finite, restricted parameter range. Using 
this approach, we first study the Hubbard model on the half-filled square lattice bilayer, including 
an interlayer repulsion term in addition to the local repulsion, and present the ground state phase 
diagram within the accessible parameter region. Starting from the antiferromagnetically ordered 
state in the absence of interlayer repulsion, the interlayer interactions are found to destabilize the 
antiferromagnetic order, leading to a band insulator state. Moreover, for sufficiently strong interlayer 
tunneling, we also observe the emergence of a direct dimer product state of mixed D-Mott and S-Mott 
character along the equal coupling line. We discuss the stability range of this state within strong¬ 
coupling perturbation theory. Furthermore, we consider the Hubbard model on the honeycomb 
lattice with next-nearest-neighbor interactions. Such an interaction is found to enhance both charge 
density and spin-current correlations within the semimetallic region. However, inside the accessible 
parameter region, they do not stabilize long-ranged charge density wave order nor a quantum spin 
Hall state, and the only insulating state that we observe exhibits long-range antiferromagnetism. 


I. INTRODUCTION 

In the field of strongly correlated electrons, the two- 
dimensional single-band Hubbard model is arguably the 
most studied model, particularly in the context of the 
metal-insulator transition and high-Tc superconductiv¬ 
ity. Even though the long-ranged Coulomb repulsion is 
considered as screened to a local interaction among elec¬ 
trons, the Hubbard model supports a variety of phases, 
including (anti-) ferromagnetism and (to a less certain 
degree) also unconventional superconductivity. Adding 
nonlocal repulsive interactions opens the large chest of 
extended Hubbard models with even richer phase dia¬ 
grams. However, in comparison to the standard Hub¬ 
bard model, such extended models have been studied far 
less by numerical means, in particular beyond one spa¬ 
tial dimension. Numerical methods that have been used 
in these instances include exact diagonalizatiorpJ^l, vari¬ 
ational cluster method^, two-particle self-consistent ap- 
proache^^*^ extensions of dynamical mean field theory 
(DMFT)l^l®or variational Monte CarlcP. 

With respect to quantum Monte Carlo (QM C) simula¬ 
tions, such studies have been rather rar c!^°fi^ l This is in 
particular due to the fermionic sign problem, which re¬ 
stricts their applicability to only small finite systems and 
high temperatures. In situations, where the sign prob¬ 
lem can be eliminated, determinantal QMC methods of¬ 
ten provide rather reliable information about emerging 
phases and the nature of the corresponding phase transi¬ 
tions. However, a direct decoupling of extended density- 
density interactions for each spin-flavor channel, such as 
employed in Refs. [10] and [TTJ leads to a sign problem even 
at half filling on a bipartite lattice, while the sign problem 
is avoided when only local interactions are considered. 
Recently, in the context of graphene, where interactions 
on the honeycomb lattice are expected to be truly long- 


ranged and thus described by a Coulomb mode l, more 
suitable QMC methods have been employed There, 
extended interactions are treated by coupling the elec¬ 
tronic density to a (continuous) scalar field, over which 
a Monte Carlo sampling is performed. Concerning the 
absence of the sign problem, this method requires - in 
addition to the usual particle-hole symmetry - the in¬ 
teraction strength to decrease sufficiently rapidly with 
distance, in a sense that will be specified further below. 

Here, we consider an alternative approach to include 
extended density-density interactions in determinantal 
QMC simulations that suffers from similar constraints 
as the aforementioned method, but which may be more 
simple to implement for models with short-ranged nonlo¬ 
cal interactions only. After introducing the approach in 
Sec. |TT] we apply our implementation to extended Hub¬ 
bard models on two different lattice geometries, both of 
which have been discussed to some extent in the recent 
literature: 

As a first application, we consider in Sec. [ITT] the Hub¬ 
bard model on the square lattice bilayer with interlayer 
repulsion at half-filling. The phase diagram of this model 
with only local repulsions has been revisited in a previ¬ 
ous publicatiorfi^, with the conclusion that the interact¬ 
ing system is either an antiferromagnetic insulator or a 
band insulator, primarily decided by the ratio of inter¬ 
layer to intralayer hopping. The addition of interlayer 
interactions was discussed previously in connection with 
exciton condensatiorP^, and has been studied with de¬ 
terminantal QMCfHl and exact diagonalizatiorP recently, 
with a focus on the impact on interlayer pairing. A recent 
cellular DMFT studj® considered attractive interactions 
and found signs of charge-ordered and superfiuid phases. 

In Sec. jlllj we investigate the effects of an interlayer 
repulsion on the phases of the half-filled system, and ob¬ 
serve the destruction of antiferromagnetic order within 
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the accessible parameter region. We then elaborate on 
the natnre of the band insulator phase and exhibit the ex¬ 
istence of an exact dimer product state inside the strong 
coupling region. We rationalise these findings based on 
large-coupling effective theorie^^^. In a related Hubbard 
model with both interlayer density interactions and in¬ 
terlayer Heisenberg exchange, a previous QMC studjff^ 
revealed the existence of a current-carrying phase. But 
also therein, the sign problem confined the explorable pa¬ 
rameter space, and the interaction range that we consider 
here has not be explored in that previous study. 

In Sec. EYi we turn our attention to the Hubbard 
model on the honeycomb lattice with next-nearest neigh¬ 
bor (NNN) interactions. Such a model has been proposed 
in a mean-field studj®! as an example for an interaction- 
driven topological Mott insulator state. According to 
that analysis, the extended interactions induce an effec¬ 
tive spin-orbit coupling, which then triggers a q uantum 
spin Hall (QSH) phase. Subsequent studies^^lMl concen¬ 
trated on the spinless variant of that model and the exis¬ 
tence of a corresponding quantum anomalous Hall phase, 
with most results pointing against its existence in that 
model, and instead finding a direct transition from a 
semimetal regime to a phase with charge density wave 
(CDW) order. Using our QMC approach, we perform 
large-scale simulations of the less-studied spinful model. 
We find that the accessible parameter range prevents us 
from observing transitions to new phases. However, we 
observe tendencies towards a CDW phase and an increase 
of QSH-compatible spin currents within the semimetal- 
lic phase prior to the instability towards antiferromag¬ 
netism. Final conclusions are presented in Sec. 


II. MODEL AND METHOD 

In the following, we are interested in the simulation of 
Hubbard models with additional nonlocal density-density 
interactions, as described by the Hamiltonian 

H = Ht + Hu + Hv, 

Ht = — 

(bV 

i 

Hv = Vijirii - l){nj - 1 ), ( 1 ) 

where rii = rii^ riii is the local density operator, U 
denotes local repulsion and the matrix elements Vij de¬ 
termine the density interaction between lattice sites i and 
j. Furthermore, denotes the free kinetic part of the 
Hamiltonian with hopping amplitudes Uj. Such models 
have been recently studied with QMC methods in the 
contex t of l ong-ranged Coulomb models on honeycomb 
lattice^i2*i^, by coupling the electron density to an aux¬ 
iliary continuous scalar field. The sign problem is ab¬ 


sent for a positive definite interaction matrix Vij (with 
Vii = U). In this section we employ an alternative ap¬ 
proach, which is particularly suited to readily include 
short-ranged interactions in auxiliary-field QMC imple¬ 
mentations that already support simulations of standard 
Hubbard models. The procedure is explained here for 
the projective auxiliary-field QMC methocP^ which al¬ 
lows for the calculation of ground state properties and 
is based on the imaginary-time evolution of a trial wave 
function |Tt) to the ground state. 


(O) 
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at which a large class of observables O can then be 
measured. The imaginary-time propagation is split up 
into M elementary increments of size At = 0/M via a 
Trotter-Suzuki decomposition. 
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The propagation by the interacting parts Hjj and Hy 
will be rendered feasible via two successive Hubbard- 
Stratonovich (HS) decouplings, 

e-5^' = — / (4) 

J — OO 


at the cost of two new auxiliary real bosonic fields 
(j)u(i,T) and over which the Monte Carlo sam¬ 

pling is eventually performed. A HS transformation re¬ 
quires the interaction term to be present in quadratic 
form A^, which in general can be achieved in different 
ways. The particular choice of HS transformations is 
crucial for the nature of the sign problerrPH and will be 
explained in the following. 

In a first step, the extended density interaction is 
rewritten in a form that couples the auxiliary field cjy 
to the particle density on the corresponding bond, 

(m - l)(nj - 1 ) =^(ni +nj - 2f - - 1 )^ 

(5) 

After identifying A? = ArVijirii + Uj — 2)^ and perform¬ 
ing the HS transformation, the propagation by Hy at 
time slice in Eq. (|^ will be replaced by 

g-Arffv (xY\ f d4>y{ij,T) e“5'^v(bT) 

i,3 

X gj -\/ArVij (Ui+Uj-2)<pv ^ 


Beside the desired quadratic term, other terms appear in 
Eq. ([^ that contribute to the local Hubbard interaction 
and hence will be absorbed into Hu. For general Vij, 
this will result in new renormalized values L/. However, 
let us assume in the following, that the extended interac¬ 
tion strength is constant and the interacting bonds form 
a lattice with connectivity zy. In this case the local in¬ 
teraction will be renormalized to U = U — zyV. 
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Since the local interaction term is already of quadratic 
form, we can straightforwardly perform a second HS 
transformation and couple the local density to the field 
</'£/{*, t), 

^-ArHu (xYlJd<j)u{i, t) 

i 

(7) 


As it turns out, the sign of the renormalized U deter¬ 
mines whether a sign problem occurs or not. More pre¬ 
cisely, sign problem-free calculations require a nonnega¬ 
tive value of U. The proof is analogous to the one given in 
Ref. Due to the SU{2) symmetry of Eq. Q, the con¬ 
figuration weight can be factored into spin components, 
WaM. {(t>v)] = n,T=td 


W,[{cl)u},{M] =Tr 


M 
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( 8 ) 


where the Hamiltonian Ht = J2ija ^la[^T]ijCja is used 
for the generation of the trial wave function, with 
Ht I'I't) = Et The sign problem is absent, if one 

can find a canonical transformation that connects both 
weights as Wi = W^, ensuring a positive total weight 
W = > 0. The SU{2) symmetry to¬ 

gether with a particle-hole symmetry allows for a canon¬ 
ical transformation —?► (—where the factor (—1) 
only contributes on one of the two sublattices. Its appli¬ 
cation to W-f and a subsequent complex conjugation will 
- except for the spin flip - not affect the trial wave Hamil¬ 
tonian and the kinetic contribution, as long as both de¬ 
scribe bipartite systems. The particle density will trans¬ 
form as ^ 1 — riii, introducing a minus sign in the 
corresponding exponentials. This minus sign will only 
get canceled by the complex conjugation, if both terms 
are purely imaginary, which is indeed the case for 

U = U-zvV>0. (9) 

In the case of nearest-neighbor (NN) repulsion only, it 
turns out that this constraint is identical to the positive 
(semi)definiteness of the interaction matrix Vij, which is 
required for the approach employed in Refs. [12] and [131 
This is no coincidence, and in fact, our approach can be 
obtained within this framework. While in the approach 
of Refs. [T^ and the scalar field couples to the largest 
modes of the global interaction matrix dynamically dur¬ 
ing the simulation, we carried out the diagonalization lo¬ 
cally, and thereby couple to these largest eigenvectors in 
terms of local, bond-centered auxiliary fields. For further 
extended interactions, such as the case of NNN repulsion 
on the honeycomb lattice considered in Sec. EYl we can 
in fact arrive at a less restricted interaction range than 
implied by Eq. by employing an appropriate local de¬ 
coupling of the interaction terms, as discussed in detail 
in Sec. |IV| Finally, we note that in practice, we employ 
a discrete version of the HS auxiliary field decoupling 


instead of the continuous field considered in the above 
derivation, e.g., the SU(2) symmetric HS decoupling of 

Ref.Hi 


III. SQUARE LATTICE BILAYER WITH 
INTERLAYER REPULSION 

In this section, we examine the half-filled square lattice 
bilayer with finite interlayer repulsion, as described by 
the Hamiltonian 

H = -tY 4xaC3X,y - + 42aC^la) 

{ij)<jX icr 

+ ^ “ ^){ni2 - 1), 

i\ i 

( 10 ) 

where A S {1,2} denotes the layer index, t± the interlayer 
hopping amplitude and V± the interlayer density-density 
repulsion. The interaction connectivity of the nonlocal 
interaction for this model is = 1, hence from the 
condition in Eq. ([^ it follows that we can simulate this 
model in the region V± < U without a sign problem. In 
contrast, in a recent QMC study of the same in model in 
Ref. [TTJ where a separate, direct HS decoupling of the in¬ 
terlayer repulsion for each spin-spin sector was employed, 
a sign problem arises for any finite values of Ul > 0 al¬ 
ready for the half-filled system. 


A. Phase diagram 

In the limiting case of vanishing in terlay er interactions, 
V± = 0, there have been suggest ion^^SM] for a paramag¬ 
netic metallic phase persisting at small values of U/t in 
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FIG. 1. (Color online) Phase diagram of the half-filled Hub¬ 
bard bilayer in the absence of interlayer interactions (Vx = 0), 
with estimates of the phase boundary between the AF insu¬ 
lator and the band insulator as obtained from QMC and fRG 
calculations. The inset illustrates the square lattice bilayer 
geometry. Figure reproduced from Ref. 1141 


the regio n of low tx/t < 4. However, more recent anal- 
yge^HEU aj.gu 0 for a direct onset of an antiferromagnetic 
(AF) insulator at any finite interaction strength in the 
low region. Increasing the interlayer hopping t^/t 
however favors singlet formation on the interlayer bonds 
and eventually drives the system into a band insulator 
phase. The Vx = 0 phase diagram as obtained from 
a combined QMC and functional renormalization group 
(fRG) analysi^Slis shown in Fig. [^for reference, while in 
the following we will study the effects of finite Vx > 0 on 
the phase diagram. 

For the calculations at finite Vx, we employed the pro¬ 
jective QMC method along with the extensions intro¬ 
duced in Sec. |TT] taking the ground state of the non¬ 
interacting model as the trial wave function. We used 
a projection length 0 = 60/t and a discretization step 
Ar = 0.1/t, chosen sufficiently low such that correspond¬ 
ing systematic errors are below the statistical uncertain¬ 
ties. The available computational resources allowed us 
to simulate systems with a maximum linear length of 
L = 24, corresponding to fV = = 1152 lattice sites. 

We restrict our analysis on the intermediate-to-strong 
coupling regime in the following, since QMC studies of 
the weak-coupling region {U < 4t) suffer not only from 
exponentially small AF order parameters in that region, 
but also from large finite-size effects due to the Fermi sur¬ 
face structure of the noninteracting tight-binding model 
on the bilayer square lattic^i^. 

Besides the AF and band insulator phases, a charge 
density wave (CDW) state with alternating electron den¬ 
sities on the two sublattices may be expected to stabilize 
for sufficiently high values V^jU. The relevant observ¬ 
ables for the ordererd states are thus the structure factors 
for antiferromagnetic and staggered CDW order, respec- 
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FIG. 2. (Golor online) Eflect of interlayer repulsion on the 
AF (left) and GDW (right) structure factors for different in¬ 
terlayer hoppings tx/t. The data shown is obtained from 
simulations of a L = 12 system at U = 8t. 



lively, given by 

^ E [(S.A • 5,a) - (5.a • Sf^)] , 

-S'cDW = E [{{nixTijx)) - {{nixrijx))] , 

( 11 ) 

where Six is the spin vector operator, Q = (tt, tt) is the 
ordering vector and A denotes the opposite layer of A, 
i.e., 1 = 2, 2 = 1. The cumulant is defined as ((OP)) = 
(OP) — (O) (P) for two arbitrary operators O, P. 

In Fig. both these structure factors are shown as a 
function of the interlayer repulsion for different interlayer 
hopping strenghts, as obtained for a system with L — 12, 
U = 8t. Increasing the interlayer repulsion can be seen to 
strongly suppress the AF tendencies, while it only weakly 
augments the CDW signal. For t± = t, a doubling of the 
CDW structure factor from Vx = 7.9t to Vx = 8 t can 
be observed, hinting at a possible first-order transition 
at this point. However, a careful finite-size extrapolation 
to the thermodynamic limit of the corresponding order 
parameters rUg = a/ Saf/N and rric = ^Scdw/N sup¬ 
ports a finite AF order parameter (Fig. and thus the 
system remains in the AF phase also at this point. In 
the opposite case of large interlayer hopping t± = 3t, the 
system resides within the band insulator phase already 
at zero interlayer interaction (cf. also Fig. [^, which is 
found to be stable for all accessible values of Vx. Only 
for the intermediate case of t± = 2t, do we observe a 
phase transition between the AF and band insulator in¬ 
duced by Vx. Performing a finite-size scaling analysis of 
the AF structure factor, with the scaling ansatz 

Saf/N = L-"-^:Fs{gLi), 5 = (12) 

Fx,c 

and using the critical exponents of the three-dimensional 
Heisenberg 0(3) universality clas^^, one can indeed lo¬ 
cate the transition point at Vx,c = 2.24(3)t. The finite- 
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FIG. 3. (Color online) Finite-size behavior of spin (disks) and 
charge (triangles) order parameters along with linear extrap¬ 
olations in 1 /L to the thermodynamic limit after the observed 
jump of the CDW structure factor at Vx = C/ = 8t, t± = t. 



FIG. 4. (Color online) Finite-size scaling analysis of the phase 
transition at [/ = 8t, t± = 2t. Left panel: Crossing plot for 
determining the critical V±. Right panel: Data collapse from 
using the critical exponents of the three-dimensional Heisen¬ 
berg 0(3) universality class. 


size data exhibit a crossing point at this critical value 
of V±, and furthermore produce a good data collapse in 
the critical region, as shown in the left and right panel of 
Fig.i respectively. Repeating the above scaling analysis 
for other values of the local interaction U results in the 
phase diagram shown in Fig. 

Within the accessible region of V± < U, we observe 
only an AF phase and a band insulator phase with a 
phase boundary that moves towards smaller interlayer 
hopping strength upon increasing U/t. This qualitative 
behavior presumably extends also into the weak coupling 
regime. 


B. Direct dimer product state for large V± = U 

While analysing the equal coupling case, Vj_ = U, we 
observed a rather peculiar behavior in the QMC dynam- 


FIG. 5. (Color online) Phase boundaries between the AF and 
the band insulator phase for different values of the local inter¬ 
action strength (lines are guides to the eye). For V± < U, the 
transition is described by the three-dimensional Heisenberg 
0(3) universality class. At V± = U, the onset of magnetiza¬ 
tion is estimated by finite-size extrapolations. 



QMC bin index 


FIG. 6. (Color online) QMC dynamics of the interlayer spin- 
spin correlations at equal coupling U = V± = 8t for ti_ = 
1.2t, 1.3t. The dashed horizontal line is set at a value of the 
ordinate of -0.375. For values t±_ > 1.3t, the ground state is 
given by the DDPS defined in Eq. (151. 


ics. We find that above a critical value of the interlayer 
hopping tx,c and after only a short thermalization period, 
all observables appear to be fixed to certain values, and 
exhibit no more fluctuation in the further course of the 
simulation. Depending on the type of observable, the ac¬ 
tual values are even independent of the choice of Ujt and 
tx/i- For example, the interlayer spin-spin and density- 
density correlations take on values Sn ■ Sa = —0.375 and 
niini 2 = 0.5, respectively. Most remarkably, all correla¬ 
tions between sites that do not share an interlayer bond 
vanish. This observation strongly suggests that in these 
cases the ground state forms a direct dimer product state 
(DDPS), located on the interlayer bonds. 

The precise nature of the DDPS can indeed be accessed 
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by an effective theory in the strong-coupling regime. 
Such a description has been employed in Ref. [T7| for the 
case of the extended Hubbard model on the two-leg lad¬ 
der, and can be readily generalized to the two dimen¬ 
sional case. In the following, the derivation will only be 
sketched and for a more detailed discussion we refer to 
Ref. [T3 

In the limit of vanishing hoppings t,t± —>■ 0, the wave 
function can be written as a product of dimer states resid¬ 
ing on the interlayer bonds. In particular for comparable 
coupling strengths, V± « U, one then finds that only 
two such states contribute to the ground state manifold, 
which are the so-called D-Mott and S-Mott states, 


( 13 ) 

While the D-Mott state describes a spin singlet, the 
S-Mott is a symmetric state with strong charge fluc¬ 
tuation. These states are now identified as the eigen¬ 
states of a pseudospin operator Sf = niini 2 — with 
Sf\D)^ = \ \D)^ and Sf\S), = -i |5),. Treating fi¬ 
nite hopping terms in second-order perturbation theory 
and subsequently performing a sublattice spin rotation, 
—>■ yields an effective two-dimensional 

XXZ spin model with a staggered transverse field and 
a uniform longitudinal field h^, 

H = + SfS^ + AS^S]) 

(ij) 

-J2{{-iyhJ: + h,s^), ( 14 ) 

i 

Q J 2 C 

with effective parameters J = ^, A = = U — V± 

and = 4t_L. 

Here, we are interested in the regime of equal coupling, 
Vj_ = U, and strong interlayer hopping t±, which trans¬ 
lates in the effective XXZ model to a vanishing transver¬ 
sal field hx = 0 and a large longitudinal field hz, re¬ 
spectively. In this limit, the XXZ spin model has been 
studied thoroughljEHMl _ in some of these references in 
an equivalent formulation as a hard-core boson model - 
and the phase diagram is well understood. For zero or 
weak longitudinal field hz , the model is in an AF ground 
state, with the staggered magnetization aligned in the 
z-direction due to the easy-axis anisotropy. Increasing 
the field strength eventually triggers a first-order spin- 
flop transition, which then saturates into a fully polar- 
ized state at hz = ^J, corresponding to t± = 4|j in the 
original bilayer Hubbard model. 

In terms of the original Hubbard model, the pseudospin 
state fully polarized in the S''^-direction corresponds to a 
product of the local eigenvectors \Sf = -I- 5 ) = {\D)^ + 
\S)y)/'/2, i.e., the bilayer Hubbard model ground state 


\D),= 




[t\ 

A 1 


tAJ 


\S),= 


V2 






with 




| 0 ). 


is given by the DDPS 




(15) 


Indeed, such a state is consistent with the observed values 
of Sii ■ Si 2 and nuna, and its energy per site of E/N = 
(VI — 2t±)l2 is reproduced by the simulations as well. 

We can now exploit the particular behavior of this 
state in the QMC simulation and scan t±/t in steps of 
At± = O.lt for various values ofV± = U to determine the 
critical hopping strength tj_^c throughout the band insu¬ 
lating phase. These calculations were performed for lin¬ 
ear size L = 8 and occasionally compared with L = 16, 24 
results, where no noticeable finite-size effects on the criti¬ 
cal hopping were observed. The resulting stability line is 
shown in Fig. and separates the DDPS from ground 
states with finite entanglement between the interlayer 
dimers. Within our accuracy, the stability line of the 
DDPS falls within the estimated transition region out of 
the AF phase (cf. Fig. |^. However, we are not able to 
discern, whether indeed a direct transition from the AF 
phase to the DDPS takes place, or whether there exists a 
small, intermediate band-insulator region that separates 
the DDPS from the AF phase. We noticed however, that 
data taken at V± = U for values of t± close to the transi¬ 
tion region within the AF phase does not lead to a data 
collapse such as those discussed in the previous section. 
This may be considered as an indication, that the tran¬ 
sition out of the AF phase at V± = U is of a different 
character than the generic transition out of the AF phase 
into the band insulating regime. However, within our 
simulations, we also did not detect clear indications for 
a (first-order) direct transition from the AF phase to the 
exact DDPS, which one may expect from the fact that 
within the DDPS the correlations vanish exactly for all 
sites not belonging to the same interlayer dimer. The ac¬ 
tual behavior within the transition region to the DDPS 
at VA = thus remains an open question for possible 
further research. 

Fig. 0 also contains the phase boundary obtained 
within the effective XXZ theory, which predicts the 
DDPS to disappear at smaller values of One needs 
to note however, that the AF state is not contained in 
the effective theory’s state manifold and thus we do not 
expect it to describe the transition correctly. 

As a general side note, we remark that it can be shown 
that the DDPS is an exact eigenstate of the Hubbard 
model for any lattice geometry constructed from dimers, 
as long as the intradimer repulsion strength coincides 
with the local repulsion strength. However, we are not 
aware of an obvious criterion that determines whether 
this eigenstate will also be the ground state. 


C. Strong interlayer repulsion 

Tuning the interlayer repulsion V(l beyond the value of 
the onsite interaction U corresponds in the effective XXZ 

















7 



0.0 -^^^^^^ 

0 2 4 6 8 10 12 

ujt = V±lt 


with the effective model parameters given as 


J = 


2t^ 

214 - t7’ 


h = 


2t‘]_ 

V±-U 


> 0 . 


(19) 


This quantum Ising model exhibits a quantum phase 
transition for the ratio of r = h/J equal to a critical 
value of Tc = 3.044(1)P^, between a low-r AF phase and 
the large-r paramagnetic phase. In terms of the original 
Hubbard model, these phases correspond to the CDW 
and S-Mott phase, respectively. Within the perturba¬ 
tive calculation, the stability line of the CDW phase is 
obtained as 


t±/t = 




v±-u 

214 - U" 


( 20 ) 


FIG. 7. (Color online) St abili ty line for the direct dimer prod¬ 
uct ground state of Eq. (15l. The red circles represent the 
QMC results, while the dashed line shows the prediction for 
the stability line obtained from the effective spin model in 
Eq. pi. 


model in Eq. (14) to a finite value of the transverse mag¬ 
netic field hx = U — V±. Starting at hx = 0 from the di¬ 
rect product state in Eq. (dn, which represents an equal- 
weight superposition of the D-Mott and S-Mott states, a 
finite transverse field thus enhances the D-Mott (S-Mott) 
character of the ground state for hx > 0, corresponding 
to U > V± {hx < 0, corresponding to V± > U). For the 
case of dominant local interactions U, this observation 
is in accord with resul ts from previous QMC studies of 
the t/-only modeP^*^, which at large identified a 
band insulator phase that adiabatically connects to the 
dimer-singlet state in the large C-limit, which is a state 
of D-Mott character. 

While we cannot perform QMC simulations in the 
large-Vj_ region, we may nevertheless employ a pertur¬ 
bative calculation to access this regime. In particular for 
the region V± ^ U ^ t,t±, after a decoupling of the 
interaction terms in the interlayer rung states the appro¬ 
priate basis for the single rung states is set by the two 
state comprising the S-Mott basis stat^i^. These can be 
expressed upon introducing a pseudospin notation as fol¬ 
lows: 


so that the CDW state is stable below t±jt = \pr\j2 « 
1.23 in the limit 14 ^ U. For beyond this value, 
the anticipated CDW state with long-range charge order 
is replaced by a dimerized state of S-Mott character. We 
expect this state to adiabatically connect to the S-Mott 
state that emerges for I 4 ^ U, and which we considered 
above. We thus find that for sufficiently large interlayer 
hopping t±/t, both the AF and the CDW state are un¬ 
stable towards dimerization into insulating states with D- 
Mott and S-Mott character, respectively. These ground 
states become degenerate and form the DDPS at equal 
coupling, 14 = U, for sufficiently large values of t±/t, 
as confirmed by our QMC simulations in the previous 
section. 


IV. HONEYCOMB LATTICE WITH V 2 
INTERACTIONS 

Here, we consider the half-filled Hubbard model on a 
honeycomb lattice with second-nearest neighbor density 
interactions, 

^ = + 4 - 1 )^ 

{ij)a i 

+ y “ 1)K' - 1)- (21) 

«b'» 


l+)i - 




(16) 


with corresponding pseudospin operators t/ rf, that 
act as follows on the above states: 


rn±),=±|±),, Tr|±), = |T),. (17) 

By a calculation similar to the case of the extended Hub¬ 
bard model on the two-leg laddeili^, one then obtains 
within second-order perturbation theory in t, t±, an effec¬ 
tive quantum Ising model on the two-dimensional square 
lattice, 

^ = (18) 
(hi) * 


where {{ij)) denotes next-nearest neighbors (NNN). 

Since the coordination number of the V 2 interactions is 
ZV 2 = 6, a direct application of the method described in 
Sec.|TT]would enable us to access a rather small parameter 
region V 2 < Uj6 only. However, the particular symmetry 
of the NNN interactions on the honeycomb lattice allows 
for a partition into N interacting triangles, each spanned 
by three lattice sites from one sublattice (cf. the inset in 
Fig.[§. This enables a rewriting of the NNN-interaction 
terms of the involved lattice sites, 

(n* - l){nj - l)+{nj - l){nk - 1) + {nk - l){ni - 1) 

= l(ni + nj +nk- 3)^ - l(ni - 1)^ 

_i(^n,-ir-l{nk-ir, ( 22 ) 
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FIG. 8. (Color online) Phase diagram of the half-filled Hub¬ 
bard model on the honeycomb lattice with next-nearest- 
neighbor interactions. The circles indicate the position of 
the phase boundary as obtained from QMC, while the solid 
line is a guide-to-the-eye. The gray dashed line marks the 
boundary of the QMC-accessible region. The inset sketches 
one three-site triangle to whose density the auxiliary field is 
coupled. 


so that a subsequent HS transformation couples the aux¬ 
iliary field (pv to the total density of each triangle. This 
type of coupling proves very beneficial as it doubles the 
sign problem-free simulation range to V 2 < U/3 and 
moreover reduces the number of auxiliary-field locations 
from 3iV to N. Note that the eased parameter restriction 
does not permit calculations beyond the applicability of 
the methods used in Refs. im and Ea as it is identical 
with requiring positive (semi-)definiteness of the interac¬ 
tion matrix. 

We implemented the above idea into the projective 
QMC method and used a projection length of 0 = 60/f 
and a imaginary-time discretization of At = 0.1/f. The 
trial wave function is generated from the noninteracting 
system’s ground state. 

The results of our investigation are summarized in the 
phase diagram in Fig.At V 2 = 0, the system undergoes 
a quantum phase transition from a semimetal (SM) to 
an antiferromagnet (AF) at Uc = 3.843(8)P^, described 
by the chiral-Heisenberg universality class of the Gross- 
Neveu-Yukawa theory^Sl. There has been some debate 
over a possible interm ediate quantum spin liquid phas^^, 
however newer result^^HlIlllSl obtained on larger clusters 
are more consistent with a direct transition. 

For all accessible finite values V 2 <Uj3 neither a QSH 
nor a CDW phase can be detected. Thus, the SM and 
AF phases remain stable, with only their phase boundary 
being shifted to larger values {Uc = 5.00(2)t on the V 2 = 
U/3 line). In the following, we report the details of our 
calculations for which this phase diagram is obtained. 



-27r -7r 0 7r 



.5 

.4 


■3 8 

Co 

.2 


.1 

.0 


qx/a 


qx/a 


FIG. 9. (Color online) Diagonal charge structure factor S{q) 
at [/ = 3t, L = 18. The peaks emerge at the FC-points of the 
Brillouin zone. Left: V 2 = 0. Right: V 2 = t. 



U/t = 3 F 2 A V2lt 


FIG. 10. (Color online) Left: Behavior of S'af (red circles) 
and Scdw(A’) (blue triangles) as a function of U at largest 
possible V 2 = t//3 for system size L = 12. Right: Comparison 
of spin and charge structure factors in the semimetal phase 
at U = 3t and L = 12. 


A. Nature of the large-F2 phase 

Recently, Raghu et al.^ suggested that sufficiently 
strong V2“ii4teractions in this model induce a phase tran¬ 
sition from a semimetallic phase to an interacting QSH 
topological Mott insulator phase. If confirmed, this sys¬ 
tem would represent a simple instance of an interaction- 
induced topological Mott insulator. Since then, the spin¬ 
less counterpart has been studied extensivelj^^flHSSJ 
most results arguing against the existence of a stable 
topological phase. Instead, at large V 2 the system ap¬ 
parently enters a charge-modulated phase, which can be 
understood as a frustrated CDW with an ordering vector 
at the Dirac point K. In the following, we will look in 
the spinful model for both, charge order and topological 
order. Charge order can be detected by monitoring the 
corresponding structure factors 

ScDwi^) = ^ [{{n,csnys)) ± , 

x,y 

s^A,B 

(23) 

where x, y run over all unit cell positions, s € {A,B} is 
the sublattice index and s denotes the opposite sublattice 
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to s. The cumulant is defined as {{OP)) = {0P) — {0){P), 
for arbitrary operators O, P. As it turns out, the effect of 
V 2 -interactions can be most clearly seen in the diagonal 
structure factor component 5 'cdw(q) = |[‘5'cdw(9) + 
‘^CDw(9 )]> which intersublattice correlations are not 
considered. In Fig. we show the diagonal charge struc¬ 
ture factor component for a L = 18 system at U = 3t, 
comparing the V 2 -free case with the system at the largest 
accessible interaction V 2 = t. While the former appears 
devoid of any features, at finite values of V 2 a honeycomb¬ 
like structure emerges with maxima at wave vector K, 
signaling a tendency towards a charge density wave with 
such an ordering vector. We thus identify Scdw{K)/N 
as an appropriate order parameter and in Fig. we 
study its behavior inside the explorable parameter regime 
alongside the antiferromagnetic structure factor, 

5] [{Sa.S-Sys)-{S^S-Sys)]. (24) 

x,y 

s=A,B 

Being most promising for the existence of a CDW phase, 
we immediately focus on the boundary of the sign- 
problem-free region and perform a scan along the V 2 = 
U/3 line. One observes that around U ~ the AF 
structure factor begins to grow sizably, while the CDW 
decays correspondingly (note the logarithmic scale), in¬ 
dicating the onset of the AF phase. In the semimetal 
region {U < At), where Scdw{K) is comparatively large, 
we can study the influence of V 2 at fixed U = 3t and find 
that increasing V 2 noticeably weakens AF and strength¬ 
ens CDW ordering tendencies. However, the AF tenden¬ 
cies remain dominant throughout and a potential phase 
transition requires larger values of V 2 , which we cannot 
access here. The finite-size analysis at U = 3t, V 2 = t 
in Fig. [TT] confirms this picture and reveals the algebraic 
behavior of all considered order parameters, and thus in¬ 
dicates that the semimetallic phase is still present here. 

Identifying a possible topological Mott insulator phase 
proves more difficult than a CDW. There exist several ap¬ 
proaches to detect the associated topological phase tran¬ 
sition in QMC simulations, involving the meas urement of 
Chern number^^, entanglement spectreP^^^ or strange 
correlators^. However, from our previous analysis we do 
not expect such a phase transition to occur in the acces¬ 
sible region, and thus these methods will be of little help 
here. 

The QSH phase on the honeycomb lattice has been re¬ 
alized by Kane and MelS^ by adding a NNN term to the 
tight-binding description arising from spin-orbit interac¬ 
tion. 



1/L 


FIG. 11. (Color online) Finite-size scaling of AF (red squares) 
and CDW (blue circles) structure factors and the maximum 
distance spin current correlation (black triangles) at U = 3t, 
V 2 = t in a double logarithmic view. An algebraic function 
/(l/L) = aL~^ can be fitted well to the data, implying the 
persistence of a semimetallic phase at this point. 



FIG. 12. (Color online) Left: Spin current-current correlation 
pattern within one sublattice for a L = 12 system at 17 = 3t, 
V 2 = t. The arrow thickness is proportional to the magnitude 
and blue arrows denote currents flowing along the arrow di¬ 
rection. The black arrow marks the reference current. Right: 
Behavior of maximum distance spin current correlations as a 
function of the extended interaction V 2 for a L = 6 system at 
U = 3t. 


ping path bends to the right or the left. This term essen¬ 
tially generates circular spin currents inside each triangle 
of a given sublattice, and stabilizes a topological phase. 
If an interaction-induced QSH phase is realized in the ex¬ 
tended Hubbard model, its onset should be announced by 
an increase of such spin currents. We therefore monitor 
the spin current correlations {js{k,l)ja{m,n)) between 
different NNN bonds, where the spin current operator, 

js{k,l) = -i^(-l)'"(c|,^q^ - c\^Cka), (26) 


= Ht +*Aso ^ (25) 

where Aso is a measure for the spin-orbit coupling 
strength, (—1)*^ introduces a minus sign for spin-down 
electrons, and Vij = ±1, depending on whether the hop- 


is measuring the spin current flowing from site k to I. 

In Fig. [T^ , the correlations within one sublattice are 
plotted for a L = 12 system as close as possible to the 
potential QSH phase at U = 3t, V 2 = t. As in the Kane- 
Mele model, the spin currents can be seen to flow circu¬ 
larly inside each triangle, however with the correlation 
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magnitude decaying strongly with distance. To quan¬ 
tify this decay, we consider the spin current correlation 
between NNN bonds at the maximum distance on the fi¬ 
nite lattice, jTmax- Similar to the structure factors, jTmax 
shows an algebraic finite-size behavior (Fig. Ill, again 
confirming the perseverance of the semimetallic phase at 
this point. In Fig. we studied jTmax Eis a function 
of V 2 U = 3t, observing a tripling of the spin cur¬ 
rent correlations in the accessible range, which is a more 
pronounced increase than shown by the CDW structure 
factor in Fig. 

In conclusion, our results are compatible with an in¬ 
stability towards a CDW as well as an interaction-driven 
QSH phase at large V 2 , but a definite statement about 
the eventually realized phase is not possible, as we cannot 
access the corresponding parameter regime. 


B. Effect of V2 on the metal-insulator transition 

While no new phases appear within the accessible re¬ 
gion, an investigation of the effect of the NNN interaction 
on the metal-insulator transition may be performed. For 
this purpose, we employed a scaling a nalysis of the stag¬ 
gered magnetization TOs = Saf /N using the standard 
finite-size scaling ansatz including leading corrections to 
scaling, 

rusiu, L) = L-^/"(l + cL-^) (27) 

with the reduced interaction u = and the scaling 

function for the magnetization, Tm{u). In the following, 
we will assume a direct transition from a semimetal to an 
antiferromagnet, described by the Gross-Neveu-Yukawa 
theory, with critical exponents v = 1.005(5), f3 = 0.74(2) 
and the correction exponent to = 0.55(4), taken from 
Ref. 

If the rescaled magnetization is plotted against the re¬ 
duced interaction u, the data for different system sizes L 
will collapse at the correct critical value Uc onto a sin¬ 
gle curve, which then corresponds to the scaling func¬ 
tion Such an analysis has been perfomed in 

Fig. [^for three different values of the extended inter¬ 
action: At V 2 = 0, V 2 = t and also along the edge 
of the QMC-accessible region, V 2 = U /3. In all three 
cases a data collapse can be achieved, the corresonding 
critical point 14 ( 14 ) is shifted to higher values with in¬ 
creasing V 2 , starting from C/c(V 2 = 0) = 3.843(8)f over 
Ce(t4 =t) = 4.25(2)t to Uc{V 2 = U/3) = 5.00(2)t. The 
resulting phase diagram is shown in Fig. Concerning 
the value of the nonuniversal correction amplitude c, we 
adjusted it such that the data collapse for V 2 = 0 repro¬ 
duces the critical interaction Uc = 3.843t from Ref. [38j 
resulting at a value c = 0.89. 

These results appear consistent with a direct transition 
described by the Gross-Neveu-Yukawa theory, which is 
also supported by a closer inspection of the scaling func¬ 
tions for the three transition points: If they indeed be¬ 
long to the same universality class, it should be possible 


to merge their scaling functions by an individual rescaling 
via aTmibu), with real numbers a,b Ri 1. This has been 
done in Fig. where the V 2 = 0 scaling functions are 
kept fixed, while the other scaling functions are rescaled 
by factors a = 1.13, b — 1.036 for V 2 = t and by a = 0.96, 
b = 1.019 for V 2 = U/3, respectively. 


V. CONCLUSIONS 

We presented a method to include nonlocal density in¬ 
teraction in auxiliary field QMC simulations, which al¬ 
lows for sign problem-free calculations inside a restricted 
parameter region, and applied it to the half-filled Hub¬ 
bard model on two different lattice geometries. In the 
case of the half-filled Hubbard model on the square lat¬ 
tice bilayer with interlayer repulsion, the ground state 
phase diagram separates into an ordered region and a dis¬ 
ordered regime, depending on the ratio of the interlayer 
to intralayer hopping. The nature of the ordered phase 
is determined by the ratio of interlayer and local repul¬ 
sions and is antiferromagnetic for dominant local 

interaction and a charge density wave state for strong in¬ 
terlayer repulsion, which we identify within perturbation 
theory. However, the type and the position of the corre¬ 
sponding phase transition into the charge ordered phase 
could not be determined from the QMC simulations, as it 
lies outside the accessible parameter regime. Within the 
disordered, band insulator phase a peculiar state occurs 
at equal coupling, V± = U, for sufficiently large interlayer 
tunneling, which is a direct product of interlayer dimer 
states of mixed S-Mott and D-Mott character. It would 
be interesting to explore in more detail the transition re¬ 
gion between the antiferromagnetic and this direct dimer 
product state, which would however require larger system 
sizes than those accessible within our investigations. 

For the Hubbard model on the honeycomb lattice with 
next-nearest-neighbor repulsion V 2 , only the semimetallic 
and antiferromagnetic phases appear inside the QMC- 
accessible region V 2 < U/3. The corresponding phase 
transition becomes delayed upon increasing V 2 and is 
consistent with recent estimates of the critical exponents 
for the chiral-Heisenberg universality class of the Gross- 
Neveu-Yukawa theory, also in the presence of finite next- 
nearest-neighbor repulsion. Regarding the nature of the 
large-V 2 phase, results obtained within the semimetallic 
phase at the largest accessible value of V 2 are compati¬ 
ble with a possible instability of the semimetal towards 
charge order or a topological Mott insulator phase. A 
definite statement will require methods that overcome 
the sign problem also in the large -14 region. 

On a final note, we would like to mention other recent 
development which have resolved the sign problem 

for simulations of fermionic models with specific extended 
interactions also in case of an odd number of fermion fla¬ 
vors, such as spinless fermion models with repulsive (at¬ 
tractive) interactions between sites belonging to opposite 
(equal) sublattices. Based on the concept of the split 
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FIG. 13. (Color online) Data collapse of the semimetal-to-antiferror nag net transition, assuming a Gross-Neveu-Yukawa critical 
theory, and using the following critical exponents from Otsuka et al.“ v — 1.005, /3 = 0.74, tj = 0.55. The amplitude c = 0.89 
has been determined such that the data in the V 2 = 0 case collapses at Uc = 3.843t, in agreement with the result by Otsuka et 
aP. Left: V 2 = 0. Center: Ya = t. Right: V 2 = U/3. 



FIG. 14. (Color online) Collapse of all three scaling functions 
from Fig. |13| The used scale factors are a = 1.13, b = 1.036 
for V 2 and a = 0.96, h = 1.019 for V 2 = U/3. 


orthogonal group, these different approaches have been 


unified, providing insightful guiding principles for identi¬ 
fying sign problem-free Hamiltonian^^, which illustrates 
the dormant potential of QMC methods. 
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